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Abstract 

The hypothesis of ecological speciation states that as populations diverge in different niches, reproductive isolation evolves 
as a by-product of adaptation to these different environments. In this context, we used Nolina parviflora as a model to test if 
this species evolved via ecological speciation and to explore current and historical gene flow among its populations. Nolina 
parviflora is a montane species endemic to Mexico with its geographical distribution restricted largely to the Trans-Mexican 
Volcanic Belt. This mountain range is one of the most complex geological regions in Mexico, having undergone volcanism 
from the mid-Miocene to the present. Ecologically, the Trans-Mexican Volcanic Belt possesses different types of vegetation, 
including tropical dry forest; oak, pine, pine-oak, and pine-juniper forests; and xerophytic scrub - all of which maintain 
populations of N. parviflora. Using species distribution models, climatic analyses, spatial connectivity and morphological 
comparisons, we found significant differences in climatic and morphological variables between populations of N. parviflora 
in two distinct Trans-Mexican Volcanic Belt regions (east vs. west). This could mean that the geographically isolated 
populations diverged from one another via niche divergence, indicating ecological speciation. Spatial connectivity analysis 
revealed no connectivity between these regions under the present or last glacial maximum climate models, indicating a lack 
of gene flow between the populations of the two regions. The results imply that these populations may encompass more 
than a single species. 
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Introduction 

The ecological speciation hypothesis postulates that incipient 
species adapt to different niches, which are influenced by the 
ecological environment (e.g. resources, predators or abiotic factors) 
[1], ultimately resulting in niche differentiation and leading to 
speciation. Svensson [1] defines three critical conditions necessary 
for ecological speciation: 1) niche differences must exist with 
sufficient magnitude between species such that divergent selection 
can drive speciation, 2) niche differences must persist for a 
sufficiently long period of time to allow reproductive isolation to 
evolve, and 3) ecological differentiation must evolve prior to 
reproductive isolation such that the ecological differentiation is 
driving the process. However, even given these criteria, ecological 
speciation is only one factor in the speciation process, with other 
geographical processes often contributing to reproductive isola- 
tion, and with confounding factors associated with ecological 
differentiation also having an effect (e.g. allopatric, parapatric, and 
sympatric modes of speciation; [2,3]). 

The ecological niche concept is most purely applied when the 
geographical ranges of the species in question overlap in 
distribution [4,5]. There are two generalizable modes of speciation 
that can occur in terms of the niche, each with testable hypotheses 



in the phylogenetic context: niche conservatism and niche 
divergence [5] . Niche conservatism is the tendency of species to 
retain ancestral ecological characteristics, and many aspects of the 
fundamental niche can be conserved over long evolutionary time 
scales [6,7]. As such, closely related species are expected to have 
more similar niches than would evolve randomly [6]. Niche 
divergence occurs when new species adapt to new niches as part of 
the speciation process, such that sister species are more divergent 
in their niche space than would be expected by chance [5] . In the 
case of allopatric speciation, niche conservatism could play a role 
early on as species could retain their ancestral niche preferences in 
geographic isolation; subsequent ecological divergence between 
the two species would limit any further gene flow between them 
[1,5] and provide mechanisms for continued isolation if their 
ranges were to overlap eventually. Thus, niche conservatism and 
ecological divergence can work together, but at different stages in 
the speciation process [1]. The existence of highly conserved 
ecological niches between closely related species indicates that 
many taxa rely on the conservation of ancestral niche space and its 
effects in causing population fragmentation and allopatry, and that 
this is an important force in speciation [8,9]. The niche 
conservatism hypothesis proposes that species distribution patterns 
are governed by ancestral climatic affinities [10]. 
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Figure 1. Geographical distribution of Nolina parviflora along the Trans-Mexican Volcanic Belt. Population numbers are those given in 
Table 1. The bottom tree shows the results of the Bayesian molecular dating analysis of N. parviflora, numbers inside each clade indicate the 
populations included in the analysis, and node numbers indicate divergence time and its 95% highest posterior density (HPD) intervals (adapted from 
Ruiz-Sanchez & Specht, 2013). 
doi:1 0.1 371 /journal.pone.0098754.g001 



Recent studies have proposed new methodologies to test the 
importance of niche conservatism and niche divergence in 
speciation [5,11,12], particularly using GIS-based environmental 
datasets. Using these methods, researchers have found evidence 
for speciation by niche conservatism in temperate montane regions 
[9,12]) and in the tropics [5,6,13]), or for niche divergence along a 
climatic gradient in the tropics [14—16]. 

Variations in climate and landscape have been shown to play an 
important role in divergence at the population level [17]. Global 
climatic cycles (e.g. glacial maxima) force populations to retract or 
spread with decreasing or increasing habitat availability, and 
different species respond differently to these changes [17]. Species 
with distributions in temperate habitats and at low elevations tend 
to retract as they lose habitat during glacial maxima [18], causing 
restricted gene flow among the remaining populations. Species 
with populations distributed at higher elevations tend to expand 
their ranges, and thus gene flow is increased among populations 
with the expansion of suitable habitat [17,19,20]). Thus, landscape 
changes in response to climate change can result in the formation 
of patches of distribution scattered among optimal habitats [17]. 
The resulting gaps in distribution lead to a reduction in gene flow 
among fragments and permit genetic isolation [21]. 

Molina parviflora (Kunth) Hemsl., is an excellent model to test if 
ecological speciation occurred via niche conservatism or niche 
divergence, and to explore current and historic gene flow among 



populations. Molina parviflora (Nolinoideae, Asparagaceae) is a 
montane species endemic to Mexico, with its geographical 
distribution mainly along the Trans-Mexican Volcanic Belt 
(TMVB) in central Mexico (Fig. 1), where it inhabits tropical 
dry, oak, pine, pine-oak, pine-juniper forests and xerophytic scrub 
from 1700 to 2800 m above sea level [22]. The TMVB is a 
complex geological region with a high degree of topographic 
complexity due to its volcanic activity from the mid-Miocene to 
the present [23,24]. Recent phylogenetic and phylogeographic 
studies show that the geological complexity of the TMVB is 
correlated with the diversification and speciation of various plant 
and animal lineages [20,22,25-27]. 

Our previous phylogeographic study of N. parviflora using plastid 
DNA markers distinguished two well-supported clades separated 
geographically along the TMVB [22]. Clade 1 comprises 
populations distributed in the western region of the TMVB, while 
Clade 2 comprises the remaining populations. Divergence time 
estimates for these clades indicate a late Oligocene to Pliocene 
divergence, coinciding with two major volcanic episodes along the 
TMVB [22]. Our goal in this study is to understand how 
volcanism influenced ecological speciation in jV". parviflora and to 
test whether volcanic events during the Last Glacial Maximum 
(LGM) modified the landscape and prevented ecological connec- 
tivity among population of this species along of the TMVB. Given 
the current distribution of populations and the ancient divergence 
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Figure 2. Statistical parsimony network and map of geographical distribution of 11 cpDNA haplotypes found in Nolina parviflora. 
H1-H11: sampled haplotypes. Solid black circles: hypothetical haplotypes. Numbers in brackets indicate individuals. Pie charts represent chloroplast 
haplotypes for each sampling locality. Population numbers correspond to those in Table 1. 
doi:1 0.1 371 /journal.pone.0098754.g002 



between the two major clades, we are particularly interested in 
investigating the role of niche divergence in incipient speciation 
and population segregation by investigating the influences of 
historical climate change, volcanism, and low dispersal ability on 
current and historical patterns of gene flow in this plant species. 

Materials and Methods 

Genetic data collection 

The DNA sequences of two chloroplast spacers \psbA-tmH 
(JX178963-JX1 78986) and taL-F (JX178987-JX179010)], as 
published previously by Ruiz-Sanchez & Specht [22], were 
sampled from 22 populations for a total of 158 individuals. 
Populations extend across the full range of Nolina parviflora along 
the Trans-Mexican Volcanic Belt (Table 1, Fig. 1). 

Genetic diversity and parsimony network 

Plastid markers (psbArtrrift and frnL-F) were combined as they 
share maternal inheritance and are inherited as a unit [28]. A 
statistical parsimony haplotype network of the populations of N. 
parviflora from the TMVB was generated with TCS 1.2.1 [29] with 
a 95 % of probability connection limit and treating gaps as missing 
data. Genetic diversity (h), nucleotide diversity (n) and the number 
of segregating sites (Sj within populations was assessed using 
Arlequin 3.1 1 [30] with the default value of 10,000 permutations. 
DnaSP v5.10.1 [31] was used to interpret pairwise F ST values 
between populations. 

Ecological niche modelling 

Species distribution models (SDM; [32]) were constructed in 
order to predict the presence of suitable habitats for jV. parviflora 
during the Last Glacial Maximum (LGM) and to investigate 
whether or not range expansion and habitat connectivity are 
observed between west and east populations of the TMVB. Ten 
environmental variables with a resolution of 1 km 2 were consid- 
ered, derived from temperature and precipitation data obtained 
from WorldClim 1.4 [33]: Mean Diurnal Range, Isothermality, 
Temperature Seasonality, Max Temperature of Warmest Month, 
Min Temperature of Coldest Month, Temperature Annual 
Range, Precipitation Seasonality, Precipitation of Driest Quarter, 
Precipitation Warmest Quarter and Precipitation of Coldest 



Quarter. These environmental variables were selected based on 
their lack of correlation (pairwise r<0.9 based on all sample 
locations) [34]. Two general atmospheric circulation models 
(GCM) were used to generate past climate scenarios for the 
LGM: the Community Climate System Model (CCSM; [35]) and 
the Model for Interdisciplinary Research on Climate (MIROC; 
[36]). Both models simulate climatic conditions as they are 
calculated to have been for the LGM, with a stronger temperature 
decrease assumed in CCSM compared to MIROC [37]. The 
original GCM data were downloaded from the PMIP2 website 
(http://www.pmip2.cnrs-gif.fr/) [38]. The models were run in 
MaxEnt 3.3.2 [39]; http://www.cs.princeton.edu/~schapire/ 
maxent/). MaxEnt employs a maximum likelihood method that 
estimates a species' distribution given maximum entropy, subject 
to the constraint that the environmental variables for the predicted 
distribution must match the empirical average [39]. 

A total of 40 georeferenced records, spatially unique and 
separated from one another by least 1 km, were compiled. 
Georeferenced data was obtained from our own collections and 
from the National Biodiversity Information Network (REMIB; 
http://www.conabio.gob.mx/remib_mgles/ doctos/ remib_ing.html; 
accessed Feb 20 1 1). The MaxEnt logistic model output for a given 
species is a continuous surface of values ranging from 0 to 1, where 
high values indicate greater habitat/ climate suitability for a given 
species. MaxEnt was run with the default settings for convergence 
threshold (10 5 ) and 1,000 interactions, ensuring only one locality 
per grid cell. To evaluate model performance, a random subset of 
25% of the total unique records was set aside, and the area under 
the curve (AUC) of the receiver operating characteristic (ROC) 
was measured providing a threshold-independent measure of 
performance. Resulting distributions were projected with QUAN- 
TUM GIS 1 .8.0-Lisboa. The output of MaxEnt consists of grid 
maps with each cell having an index of suitability between 0 and 1 . 
Low values indicate that the conditions are unsuitable for the 
species to occur, whereas high values indicate that the conditions 
are suitable. 

Climate analysis 

From the 40 unique georeferenced records, we conducted two 
different analyses; (a) multivariate analysis of variance (MANOVA) 
and (b) the multivariate method introduced in McCormack et al. 
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Figure 3. Potential distribution of Nolina parviflora identified using species distribution modelling (SDM). A = prediction of suitable 
habitat in the current environment. B = prediction projected onto the past climatic layers (LGM; CCSM), and C = prediction under past climatic 
conditions (LGM; MIROC). Blue colour indicates high probability of prediction of suitable habitat and red colour indicates areas with unsuitable 
habitat. 

doi:1 0.1 371 /journal.pone.0098754.g003 



(2010[12]). A MANOVA was used to test for differences in climate 
variables among regions of the TMVB, specifically testing for 
differences between the habitats occupied by eastern vs. western 
populations. Individual observations were binned into separate 
groups corresponding to the two clades identified in our previous 
study [22]. Group 1 ( = Clade 1) contains only populations from 
the western part of the TMVB, and Group 2 (= Clade 2) contains 
all populations located in the eastern part. To examine the overall 
level of divergence in environmental space among the popula- 
tions/regions, a Principal Components Analysis (PCA) was 
conducted using extracted values of the 10 climate variables 
(WorldClirn 1.4; Hijmans et al., 2005[33]) for each unique locality 
using Quantum Gis v. 1 .8.0-Lisboa. A MANOVA of PCA scores 
(dependent variables) was then performed to test for significant 
differences among the PCA scores and groups (fixed factors). Post 
hoc tests (e.g., Tukey's HSD) were used to detect if the regions 
differed significantly from one another, and if so, for which 
variables. The PCA and MANOVA analysis were carried out in R 
[40]. 



The McCormack et al. [12] method was used to test for niche 
divergence/conservatism of the eastern vs. western populations of 
M. parviflora. Analyses were performed with QUANTUM GIS 
1 .8.0-Lisboa. For the McCormack analysis, the occurrence data 
for two clades were first binned into separate groups, and the 
values of 10 climate variables were extracted from the occurrence 
points for each of the lineages. Then, 1,000 random points were 
generated from within a minimum convex polygon drawn around 
the occurrence sites, and values of 10 climate variables of the 1,000 
random points were extracted as the background predictions from 
within the geographic ranges of each lineage. A PCA was 
conducted on these data, extracting the first four PC (niche) axes 
for further consideration. Niche divergence or conservatism was 
evaluated on each niche axis by comparing the observed difference 
between the means for each lineage on that axis to the mean 
difference in their background environments on the same axis. A 
null distribution of background divergence was created by 
recalculating the score of background divergence over 1,000 
jackknife replicates with 75% replacement. Significance for 
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Table 3. Loadings of the environmental variables for each PC axis and tests of niche divergence and conservatism. 







PCI 


PC2 


PC3 


PC4 


West vs. East of the TMVB 


BI02: Mean Diurnal Range 


0.237 


0.354 


0.438 


-0.337 


BI03: Isothermality 


-0.301 


0.371 


0.260 


-0.265 


BI04: Temperature Seasonality 


0.430 


-0.202 


0.434 


-0.240 


BI05: Max Temperature of Warmest Month 


0.442 


—0.1 78 


— 0.338 


— 0.477 


BI06: Min Temperature of Coldest Month 


0.121 


-0.675 


-0.505 


0.202 


BI07: Temperature Annual Range 


0.421 


0.118 


0.276 


-0.139 


BI015: Precipitation Seasonality 


0.400 


-0.170 


0.159 


-0.106 


BI017: Precipitation of Driest Quarter 


-0.181 


-0.500 


0.186 


-0.357 


BI018: Precipitation of Warmest Quarter 


0.295 


-0.334 


0.265 


0.174 


BI019: Precipitation of Coldest Quarter 


-0.568 


-0.266 


-0.183 


-0.112 


Percent variance explained 


43.42 


28.63 


17.78 


3.67 


Observed difference 


3.9701 » 


1.7427* 


0.6412* 


0.2398* 


Null distribution 


4.3129-4.3145 


0.8059-0.8087 


0. 2900-0.291 1 


0. 0976-0. 0986 



* Significance level, P<0.05. Bold values indicate niche divergence. 

Observed differences in climatic niche for western and eastern lineages of Nolina parviflora on each PC compared to the middle 95th percentile of a null distribution of 
the differences between their environmental backgrounds. 
doi:1 0.1 371 /journal.pone.0098754.t003 



rejecting the null was evaluated at the 95% level. All analyses were 
conducted using Stata 10 [41]. 

Spatial connectivity 

To generate spatial resistance distances to dispersal among sites 
based on habitat suitability maps (conductance grid), we imported 
each SDM (current, CCSM and MIROC) into Circuitscape 3.1 
[42]. This program calculates pairwise resistance to gene flow 
among populations based on all possible paths, not just the least 
cost path, thus better explaining the movement of genes among 
widely separated regions over many generations [43,44]. The 
input for Circuitscape is a raster data set (habitat map) in which 
each cell is assigned a conductance value corresponding to the 
probability of the study organism moving through the habitat type 
encoded by the cell. We chose a conductance grid in which higher 
cell values denote greater ease of movement, and applied a 
connection scheme that allowed gene flow among the four nearest 
cells. We used the SDM rasters as input maps to quantify pairwise 
resistance distance between N. parviflora populations across its 
distribution along the TMVB. 

To evaluate the influence of climatic (present day, CCSM and 
MIROC niche models) factors on jV. parviflora connectivity and 
genetic structure, pairwise resistance values were generated for our 
three models (present, CCSM and MIROC) given the pairwise 
plastid genetic distances (F ST /( 1 -F ST ) [45]; Table 2). A Mantel test 
[46] was conducted to detect any significant associations between 
plastid pairwise genetic distances and pairwise resistance distances 
[17]. Mantel correlation coefficients (r) were estimated using 
IBDWS 3.16 [47] with 10,000 permutations. 

Morphological comparisons 

Morphological comparisons between east and west clades were 
performed using exemplars from the 22 populations sampled 
along the TMVB (Table 1, Fig. 1). In addition to specimens 
collected specifically for this study, herbarium material was 
examined from IBUG, IEB and XAL. Vegetative and floral 
characters were measured from living material in the field and 



from our own collections and historical herbarium specimens. 
Twelve morphological characters noted to be variable among 
individuals were scored and compared among populations. 

Results 

Genetic diversity and parsimony network 

The total length of the aligned combined plastid markers (psbA- 
tmH and taL-F) was 987 bp with 18 polymorphic sites (S). 
Haplotype diversity (h) ranged from 0-0.866, reflecting the 
presence of different haplotypes within each site (Table 1). 
Nucleotide diversity (n) was low (0.0000-0.00482), indicating little 
variation between sequences from the same population (Table 1). 
Indels of 1-5 bp were considered missing data. In general, pairwise 
comparisons of Fst values were high and most significant between 
the west vs. east/easternmost and east vs. easternmost populations 
of the TMVB (Table 2). 

The statistical parsimony network recovered 1 1 different 
haplotypes among N. parviflora populations. The most frequent 
haplotype (H6) is shared by populations from Puebla, Tlaxcala 
and Veracruz, and represents 44.3% of all individuals sampled. 
This haplotype is geographically distributed in the easternmost 
part of the TMVB. The second most frequent haplotype (HI) is 
shared by populations in the states of Edo. Mexico, Hidalgo and 
Michoacan, and represents 32.2% of all individuals sampled. HI is 
geographically distributed in the eastern region of the TMVB. 
Haplotypes H2 and H3 are shared between the Jalisco and 
Zacatecas populations in the western region of the TMVB. 
Haplotypes H4 and H5 are unique to the MICH and MEX1 
populations, respectively. Haplotype H7 is shared by the PUE2, 4, 
9 and VER populations; and finally haplotypes H8 and Hll are 
unique to the PUE4 and ZAC populations (Table 1, Fig. 2). 
Haplotypes H2-3, H9-1 1 are exclusive to the western part of the 
TMVB and are not shared with any populations occurring in the 
eastern area of the TMVB. Haplotypes HI, H4 and H5 are 
distributed in the eastern region of the TMVB and are not shared 
with any of the easternmost populations. Finally H6, H7 and H8 
are only found in the easternmost part of the TMVB (Fig. 2). 



PLOS ONE | www.plosone.org 



7 



June 2014 | Volume 9 | Issue 6 | e98754 



Ecological Speciation in Nolina parviflora 




103 102 



Figure 4. Connectivity maps among populations (black circles) for Nolina parviflora. A = current environment, B = Last Glacial Maximum 
under the CCSM model, and C = Last Glacial Maximum under the MIROC model. Blue and light shades of blue indicate areas with higher current 
density; areas where connectivity is most tenuous are shown in yellow to red colours. 
doi:1 0.1 371 /journal.pone.0098754.g004 



Ecological niche modelling 

The modelled distribution (AUC = 0.991) corresponds to the 
known distribution of Nolina parviflora in the TMVB. The SDM for 
the ten current climate variables across all populations of N. 
parviflora predicted an accurate distribution of N. parviflora for the 
eastern region of the TMVB, but underpredicts the occurrence of 
populations in the western region of the TMVB (Fig. 3A). When 
the models were projected onto past climate (2 IK) layers, two 
different scenarios were retrieved (Fig. 3B, 3C). For the climate 
layers based on CCSM (Fig. 3B), a large area of suitable habitat is 
recovered in the east and easternmost region of the TMVB and 
extends to the western area of the TMVB beyond the current 
occurrence. Habitat with some suitability connects these areas. For 
the projections under MIROC (Fig. 3C) a much smaller area of 
suitable habitat is predicted and is located mostiy in the eastern 
region, with no connection between the eastern and western 
predicted areas of suitable habitat (Fig. 3). 



Climate analysis 

The PCA of the 10 climate variables found that PC 1 = 43.42%; 
PC 2 = 28.63%; PC 3 = 17.78% and PC4= 3.67% together explain 
93.5% of the variability. Three temperature and two precipitation 
variables define PC 1 (Table 3); PC2 is defined by Mean Diurnal 
Range, Isothermality, Min Temperature of Coldest Month, and 
Precipitation of Driest and Warmest Quarter; PC 3 is defined by 
Mean Diurnal Range, Temperature Seasonality, and Min 
Temperature of Coldest Month; and PC4 is defined by Max 
Temperature of Warmest Month and Precipitation of Driest 
Quarter (Table 3). The PCA/MANOVAs that were performed 
given two separate regions (West TMVB and East TMVB) 
detected a statistically significant difference between eastern and 
western populations (Wilks' lambda =0.561; -F^s = 6.82; 
P= 0.0003). This difference occurs along PCI (F= 14.34, 
P= 0.0005) and PC 2 (77=4.96, P=0.03) indicating that these 
two regions occur in areas with differences in Mean Diurnal 
Range, Isothermality, Temperature Seasonality, Max Tempera- 
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Tests of niche divergence and conservatism of the four niche 
axes indicate niche conservatism on niche axis 1(43.42%; Table 3) 
and niche divergence for axes 2, 3 and 4 (50.08%; Table 3). 

Spatial connectivity 

Climatic spatial resistance surfaces (present, CCSM and 
MIROC) showed connectivity corridors between neighbouring 
populations in the easternmost part of the TMVB (populations 4— 
6; 7-10; 11,12,13-22) for present day conditions (Fig. 4A) and 
during the last glacial maximum (CCSM; Fig. 4B, MIROC; 
Fig 4C); but no connectivity was detected between the populations 
of the eastern and western regions, and likewise no connectivity 
was maintained between the populations of the western region 
(Fig. 4) for any of the three climatic spatial resistance surfaces 
(Fig. 4A,B,C). Spatial genetic correlations showed significant 
positive correlations between both genetic Fst/(1--Pst) pairwise 
distances and spatial resistance distances among populations in the 
current or LGM models (present: r=0.27; P = 0.005; CCSM: 
) = 0.52; P= 0.001 MIROC: r=0.28; P= 0.004). 

Morphological comparisons 

Twelve vegetative and floral characters noted to vary among 
individuals within and between populations were scored (Table 4, 
Fig. 5). When clustered into separate groups, the most significant 
differences were detected between specimens of the east v. west 
clades in (1) plant height, (2) leaf length and (3) fruit size (Table 4). 

Discussion 

In this study we integrated analyses of genetic diversity, 
ecological niche modelling, climate prediction/hindcasting and 
spatial connectivity to understand how climate, landscape and 
historical volcanism have influenced ecological divergence and 
genetic isolation in N. parviflora populations distributed along the 
TMVB. We found a clear east-west pattern of haplotype 
distribution, where haplotypes H2-3, 9-11 are found exclusively 
in the western region of the TMVB and the remaining haplotypes 
are spread throughout the eastern region (Fig. 2). The eastern part 
of the TMVB is separated by two distinct subregions; the 
easternmost populations sharing haplotypes H6, 7 and 8 and the 
eastern populations sharing haplotypes HI, 4 and 5. These two 

Table 4. Vegetative and floral morphological characters scored on specimens from the east and west clades of N. parviflora along 
of the Trans-Mexican Volcanic Belt. 



Character 


East clade average 


West clade average 


Total plant height (m) 


4-5 


1.5-2 


Diameter (cm) 


10-20 


7-9 


Number of branches per plant 


1 to 6 (8) 


1 or 2 


Leaf length (cm) 


100-150 


80-90 


Leaf width (mm) 


10-20 


8-9 


Main flowering branch length (cm) 


20-30 


9-12 


Female tepal length (mm) 


2.5-4 


2.2-2.5 


Female tepal width (mm) 


1.5-2.5 


1.8-2 


Fruit length (mm) 


7-12 


6-6.3 


Fruit width (mm) 


9-14 


7-8 


Seed length (mm) 


3-4 


3.5-4 


Seed width (mm) 


3-3.5 


2-2.5 



doi:1 0.1 371 /journal.pone.0098754.t004 
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West lineage East lineage 



Figure 5. A-C west lineage. A. Plant showing two branches. B. 
Primary inflorescence branch, showing fruits and seed, scale bar = 
1 cm. C. Panoramic image of the oak forest type in the locality of Rio 
Los Patitos, Zacatecas, Mexico. D-F east lineage. D. Plant showing seven 
branches and inflorescences. E. Primary inflorescence branch, showing 
fruits and seeds, scale bar = 1 cm. F. Panoramic image of xerophytic 
scrub habitat type with Pico de Orizaba in the background, Veracruz, 
Mexico. Photos by E. Ruiz-Sanchez. 
doi:1 0.1 371 /journal.pone.0098754.g005 

ture of Warmest Month, Min Temperature of Coldest Month, 
Temperature Annual Range, Precipitation Seasonality, Precipita- 
tion of Driest, and Warmest and Coldest Quarter. 
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subregions do not share haplotypes between them (Fig. 2). The 
high values of the F sx pairwise comparisons between populations 
from the east v. west of the TMVB indicate that there is no gene 
flow between these populations, nor between populations of the 
two subregions of the eastern TMVB (Table 2). Two recent studies 
on different species of plants with populations distributed along the 
TMVB found the same pattern of east-west isolation [48,49] . 

Ecological speciation 

Ecological niche models (ENMs) can be used to identify areas of 
suitable habitat for JV parviflora (Fig. 3). ENMs can be projected 
onto historic climate maps to predict suitable areas for ancient and 
ancestral distributions, assuming similar habitat preferences were 
exhibited by historic populations (Fig. 3B, 3C). The ENM for JV. 
parviflora given the current distribution underpredicts suitable 
habitat for some of the known and sampled populations in the 
western region of the TMVB, even though these populations were 
used to build the model (Fig. 3A). This could provide evidence for 
differential habitat preferences among populations, such that 
segregated populations occupying a different climatic niche could 
be designated as separately evolving lineages if gene flow is 
completely lacking between these lineages and those that occupy 
more central climatic niche space [7]. Another hypothesis is that 
some ecological factors (e.g. different type of vegetation) may allow 
for adaptation to different climatic niche spaces, ultimately driving 
speciation [9,10]. We observed in the field that JV. parviflora 
inhabits very different types of forests, from pine-oak-juniper 
forests to tropical dry forest, to xerophytic scrub and pine-oak 
forests. The occupation of different habitats could limit dispersal 
among populations [50] even if the populations are geographically 
proximal. 

The results of the MANOVAs suggest the occurrence of niche 
divergence between the populations of the eastern Trans-Mexican 
Volcanic Belt and those of the western TMVB. In particular, 
variables associated with temperature and precipitation (PC 1 and 
PC 2) differentiate these two regions (Table 3). The populations of 
the western TMVB are found in tropical dry and oak forests, while 
those of the eastern TMVB are located primarily in pine-juniper- 
oak forests and xerophytic scrub. Tests of niche divergence vs. 
niche conservatism indicate that half of the variation (50.08%) is 
explained by niche divergence, in agreement with results of the 
SAMOVA analysis. In addition to the complete lack of shared 
haplotypes, the high F ST values for the pairwise comparisons 
between populations from the east v. west of the TMVB indicate 
that there is no gene flow between these two regions. Furthermore, 
the SAMOVA and test of niche divergence indicate that the two 
clades could be diverged via ecological speciation as they clearly 
fulfil Svensson's [1] criteria of ecological speciation. This observed 
pattern could be the result of historical divergence, where the 
genetic drift may have acted in the smaller western populations, 
allowing divergent selection to act more rapidly given the high 
degree of reproductive isolation. 

Climatic spatial resistance surfaces for the present and the Last 
Glacial Maximum showed some corridors between neighbouring 
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populations in the eastern region of the TMVB, but no 
connectivity was maintained with the populations in the western 
region of the TMVB (Fig. 4A, B, C). The spatial genetic 
correlations showed a positive and significant (present: r=0.27; 
P= 0.005; CCSM: r=0.52; P= 0.001 MIROC: r=0.28; 
P— 0.004) correlation between genetic Fst/(1-F S t) pairwise 
distances and spatial resistance distances among populations, both 
present-day and in the LGM, indicating that the populations have 
been isolated since the LGM (Ruiz-Sanchez & Specht, 20 13 [22]) 
with no gene flow between them, except for immediately adjacent 
populations (Fig. 4). 

These results contradict those reported by Wiens [9] and Velo- 
Anton et al. [17], where the LGM permitted a continuous 
distribution that could increase the gene flow among populations. 
Instead, our results are more similar to those reported by Bryson & 
Riddle [27], where tropical montane specialists persisted in 
isolation over longer periods of time, corroborated by empirical 
studies of vegetation change in central and western Mexico over 
the last 86,000 years [51,52]. Our previous study [22] revealed 
relatively old divergence times between the eastern and western 
clades (see their Fig. 3) indicating that the uplift of the TMVB at 
different times in the past may have influenced the genetic and 
ecological diversification of JV parviflora. The Neogene formation of 
the TMVB with different episodes of volcanism [23,24] appears to 
have promoted diversification via habitat fragmentation and the 
resulting isolation of populations for different plant and animal 
lineages [20,25-27,53,54]. 

Our results provide evidence for differentiation among popu- 
lations reflecting habitat preference among extant populations, 
and some of these segregated populations could be designated as 
separately evolving lineages given that gene flow appears to be 
completely lacking [7]. Their occupation of different habitats 
could further limit dispersal between distant populations [49] , even 
more than would be expected as a result of distance alone. 
Ecological divergence could represent an important stage in the 
process of lineage separation and speciation [1,55]. Ecological 
speciation could be the key process in the speciation of JV. parviflora 
along the Trans-Mexican Volcanic Belt, a conclusion supported 
by our ecological analyses and the observed morphological 
differences identified between the eastern v. western lineages. A 
complete taxonomic description of the new Molina species is 
currently in preparation (Ruiz-Sanchez et al. in prep.). 
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